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Abstract. We introduce a simple family of models for representing the dark matter in galaxies. The potential and phase space 
distribution function are all elementary, while the density is cusped. The models are all hypervirial, that is, the virial theorem 
holds locally, as well as globally. As an application to dark matter studies, we compute some of the properties of y-ray sources 
caused by neutralino self-annihilation in dark matter clumps. 
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1. Introduction 

It is worthwhile to find simple models for the distribution of 
dark matter in galaxy haloes. This is tantamount to solving the 
collisionless Boltzmann and Poisson equations for the distri- 
bution function (DF) /, potential if/, and density p of the dark 
matter particles. 

Many authors start by assuming a profile for the dark mat- 
ter density and then solving for the self-consistent potential and 
DF. This has provided some widely-used and notable mo dels 
for dark matter haloes (e.g., Jaffe 1983; Hernquist 1990). A 
drawback to this approach is that, even if the density and po- 
tential are simp l e, the DF is often unwieldy. For example, the 
DF of theH IS modeH^hmher transcendental func- 
tion, while the DF of the Hernauist ( 1990) model is composed 
of an unwieldy bunch of elementary functions. 

It can be advantageous to tackle this problem the other way 
around, by assuming a simple DF and then solving for t he self - 
consistent potential and the density. Except for Toomre ( 1982), 
this reverse approach has not been widely used. Here, we ex- 
ploit it to build a flexible family of cusped dark matter halo 
models with an elementary DF and potential (Sects. |2]&|3j. 

As a simple application, we use our models to study the sig- 
nal from indirect detection experiments (Sect.|4j. In particular, 
y-rays from dark matter annihilation may be identified by forth- 
coming atmospheric Cerenkov telescopes such as VERITAS 1 
or by satellite-borne detectors like GLAST 2 , and so it is useful 
to have definite predictions from halo models. 



* Current address: MIT Kavli Institute for Astrophysics & Space 
Research, Massachusetts Institute of Technology, 77 Massachusetts 
Avenue, Cambridge, MA 02139, USA 

1 http://veritas.sao.arizona.edu 

2 http://www-glast.stanford.edu 



2. A Simple distribution function 

2.1. An ansatz 

Let us assume that the dark halo is spherical, in which case the 
DF may depend on the binding energy E and the magnitude 
of the angular momentum L. Let us note that the genera lized 
Plummer models, recently studied bv lEvans & AnI (|2Q05) have 
very simple power-law DFs of the form 



f(E,L)oz L p - 2 E i3p+m . 



(1) 



This suggests that it may be worthwhile to look for models with 
DFs given by the sum of such components (c ,f . . iFrickelll 952l 
lToomrell982l) : 



f(E,L)=Y J CiL p >- 2 E 0p " 



l)/2 



(2) 



where C, and are all constants. Then, the density is of the 
form 



P 



where the constants D, are 



r(p,-/2)r(3p,-/2 + 3/2) 
T(2pi + 2) 



(3) 



(4) 



By integrating the DF over velocity space, it is straightforward 
to derive the radial and tangential velocity second moments 
(v 2 ) and (v^) respectively. Whereas the anisotropy parameter 
/3 = 1 — (v 2 )/(2(v 2 )) is no longer constant in contrast to mod- 
els of Ev ans & AnI (120051) . we still find the remarkably simple 
relation between the three-dimensional velocity dispersion and 
the potential; 



<v 2 > + (v 2 T ) 



(5) 
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In other words, the root mean square velocity is always one- 
half of the escape velocity at every spot, or the virial theorem 
hold s locally for a ny model described by the DF of Eq. 0- 
lEvans & Anl |2005) coined the term "hypervirial" to describe 
such stellar dynamical models for which the kinetic energy 
in each volume element (T = p(v 2 )/2) is exactly one-half of 
the magnitude of the local contribution to the potential energy 
by the same volume element (\W\ = This idea can be 

trace d back to the c lassical investigations of Plumm3 dl91 ll) 
and lEddingtonl II 1 9 1 611 . It has received additi onal im petus from 
the recent /V-body simulations of lSota et"aD(l2005|) . 



Here, we note that the differential equation (y') 2 = f(y) where 
/(y) is a polynomial of y, the degree of which is at most 
two, can be solved through elementary functions. However, the 
right-hand side of Eq. Jl 2i can be a quadratic polynomial of 
ft if p = pi or p = 2pi for all distinct /?,'s. For the simplest 
case when the sum contains only a single term, it is straightfor- 
ward to show that either choice of p leads to the solution that 
reduces to the models of lEvans & Anl J2005I) with the integra- 
tion constant being the scalelength. The only other possibility 
is that the sum contains two terms with p — p\ — 2p 2 , which is 
investigated in the following section. 



2.2. Poisson's equation 



2.3. Solutions 



While we have established that any spherical system described 
by a DF of the form of Eq. is hypervirial, we still have 
to find the corresponding density and potential by solving 
Poisson's equation (here G — 1) 



The order of Eq. O can be reduced as follows (c.f., 
IChandrasekhail IT939h . First, let us consider the substitution 
ipexp(-t/2) and t = lnr. Then, the left hand 



r -i/2, 



side of Eq. (|5Jl transforms 



1 d 



r 2 dr 



dr 



-5/2 



d 2 tp ip 
df 2 " ~ 4 



(7) 



r- 5 V» +1 . 



while the right hand side can be rewritten using 

i i 

Hence, Eq. reduces to 

Z± = t - 4jt, V Di<p 
dt 2 4 A-> 



2p,+l 



(9) 



which does not involve the independent variable explicitly, and 
so its order can be reduced by standard techniques (e.g., Ince 
1 1944 to give 

x 2 



(dtp) 

\dt 



4 



1-2^ 



2p, 



where A is a constant of integration and 
4jtD, 

Bi 



(10) 



(11) 



Pi + 1 

Using the boundary condition at infinity (ip — and d<p/dt = 0) 
0. 3 Then, after introducing a further transfer 
ft, Eq. fllOi reduces to 



we find that A 
mation of the variable, <p 

,2 



1™) Vab^-p - 

p dt) Z-i 



(12) 



3 The solution with A 4- is unphysical because Eq. 1 101 implies 
that (d In i^/d In r) diverges as rif/ 2 — > unless A = 0. Nevertheless, we 
note that it is possible to find explicit solutions with A j= if the sum 
contains a single term with p = 1/2, 1, or 2 (See AppendixP. The re- 
sulting solutions, expressible in terms of Jacobi or Weierstrass elliptic 
functions, however, exhibit oscillatory behaviour along the real axis, 
which we suspect to be a general property of the differential equation. 



For this case, the integration results in 



ft = 2k cosh 



~(t-t ) 



+2B 2 = k 



r \ p/2 I r x ~ p/2 



ro j \ro 



+2B 2 ,(13) 



(6) where to = In ro is the integration constant, and k = {Bi+B 2 ) 1 ^ 2 . 



By reinstating ft = <p p = r p ^ 2 t// p , we obtain 



— = kr; 

i/fP 



P/2 



1 + 2c I — 

ro 



P/2 



ro 



(14) 



where c = Bijk. With the normalization M«, = 1, we find that 
k = r^ 2 . So, we arrive at a two-parameter - c and p - potential 
of the form of (incorporating the dimensional constants G, M 
and ro) 



(8) iff = 



GM 



MP' 



(r p +2crfrP' 2 + rP) 
whose DF is given by 

f(E,L) = Cl L p - 2 E 3p ' 2+1/2 + C 2 L p/2 - 2 E^ /4+ " 2 . 



(15) 



(16) 



Here, the constants C, can be found from Eqs. and Jl It 
with B\ — 1 - c 2 , Z?2 = c, p\ — p and p 2 = p/2 (henceforth 
G = M = r = 1), that is, 



C 2 



l-c 2 



F(2p + 3) 



2P' 2+1 (2n) 5 ' 2 T(p/2)T(3p/2 + 3/2) 

c Tip + 3) 

2W4+1(2jt)5/2 T(p/4)T(3p/4 + 3/2)' 



In order for the DF to be non-negative everywhere, we must 
have p > and < c < 1 . In particular, if c — or c — 1 , the 
models reduce to those o flEvans & Anl (120051) . We note that the 
model with (p, c) = (a, 1) is identical to the one with {p, c) = 
(a/2,0). 

3. The Simple halo models 

Thusfar, we have obtained the gravitational potential (Eq. 1151 
corresponding to the simple DF (Ea.ll6>. Various limits are al- 
ready wel l-known. For exa mple, when (p, c) = (1,0) or (2, 1), 
this is the Hernauistl Jl990i) potent ial generated by the DF first 
found bv iBaes & Deionghd J2002I) When (p,c) = (2,0) or 
(4, 1), this is the isotropic IPlummen i 1 9 1 1 ) model. Bearing in 
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Fig. 1. The density profile of the generalized hypervirial mod- 
els. The Plummer and Hernquist models are included in the 
generalized hypervirial family as the cases (p, c) = (2, 0) [or 
(4, 1)] and (p,c) = (1,0) [or (2, 1)] respectively. Note that the 
model with {p, c) = (1/2, 0) is the same as (p, c) — (1, 1). 

mind the property in Eq. l|5}, we refer to this family as the gen- 
eralized hypervirial models. 

The density generated by the potential of Eq. J15i is 



(l-c2)(p+l)(4ji)- 



c(p+2)(8ji)- 



(17) 



r 2_p (l +2cr p l 2 +r p ) 2+l l p r 2 ~ p / 2 (l+2cr p / 2 +r p ) 1+1 t p ' 

which can be found from Poisson's equation. Fig. [0 shows 
some typical density profiles. Provided that c + 0, the sec- 
ond term in Eq. dl7> is dominant when r — > [p ~ r~ ( - 2 ~ p l 2) \ 
and r — > oo [p ~ r -( 3 +/V 2 )]. Jf p < 2, the density is monoton- 
ically decreasing outwards-radially regardless of the value of 
c. If 2 < p < 4, there may be a region of increasing density 
depending on the value of c, although the model still exhibits a 
cuspy centre. The p = 4 model is cored whereas there is a hole 
at the centre if p > 4. 

Notice that the density profile - although composed of en- 
tirely elementary functions - is a bit more complicated than 
either the potential or the DF. We argue that this is the right 
way around as most applications will use the potential (for ex- 
ample, for integrating the orbits in numerical simulations) or 
the DF (for example, for calculating the flux of dark matter 
particles on a detector). It is much more useful to have models 
with simple potentials and DFs than those with simple density 
profile. 

Evidence from N-body simulations suggests that the den- 
sity profile of the dark halo follows a simple functional 
form . One of the most common l y cite d examples is that 
of iNavarro. Frenk. & White! 1 19951 Il996l henceforce NFW), 
which is basically a double power-law characterized by r _1 
cusp at the centre and r~ 3 fall-off at large radii. Since every 



Fig. 2. The circular velocity curves of the generalized hyper- 
virial models as a function of r/rj/2, where n/2 is the half-mass 
radius. Here, v 2 , = GM/r^. The line types are as in Fig. 1. 

member of the generalized hypervirial models has a finite mass, 
none of them can reproduce the r -3 density fall-off - which 
implies an infinite mass - in the outer region. Regarding the 
behaviour in the inner region, however, many members of the 
family indeed e xhibit a r~ l -like c usp, including the well-known 
example of the lHernquistllll990l) model. In fact, the additional 
freedom afforded by the parameter c admits more flexibility in 
the behaviour around the scalelength. For example, we find that 
the model with (p, c) = (2,3/4) provid es a better fit to the NFW 
profile within a scalelength than the iHernquisl Jl990h model. 
Actually, if we allow a slight deviation of the cusp slope, there 
exists a trade-off between varying p and c which produces very 
similar behaviour of the density profiles in the inner parts. 
The circular speed and the cumulative mass can be found 

as 



2 _ _ 
Vc r dr 



r P/2 



+ r p 



M r 



(1 + 2cr p l 2 + r p ) l l p 

1 + cr' p/2 
(1 + 2cr- p > 2 + r- p y/ p+ 



+i ' 



(18) 



(19) 



Here, the total mass is finite and therefore the circular speed 
falls off as Keplerian at large radii (v c ~ r~'^ 2 ). Fig. |2 shows 
plots of the circular velocity as a function of r/ri/2, where ri/2 
is the half-mass radius, (i.e., M n/2 = 1/2). The models with 
inner density slopes p cc r~ l or p oc r ~ 3 / 2 have ro t ation curves 
similar to that of the NFW profile or Moo re et al.l (llJ98) mod- 
els, and so are flatfish over a wide range of radii. 
The velocity dispersions are 



1 + c 2 + eg 



2(1 + c 2 + eg) + p(2 + eg) 



(20) 
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Fig. 3. Contour plots of the (projected) number density of y-ray sources in equatorial coordinates. The left hand panel refers 
to a model with p oc r _1 at the centre (p,c) = (2, 1/2), while the right-hand panel to a model with p oc r~ 7 ^ 4 at the centre 
(p, c) = (1/2, 1/2). The grey-scale gives the relative number density, with white representing the highest values and black the 
lowest. The number density changes by a factor of 10 on moving from one contour to the next. The Galactic centre and the 
Galactic plane are marked. 



(v 2 T ) 



2 + eg 



2 2(1 + c 2 + eg) + p(2 + eg) 



(21) 



where g = r p ^ 2 + r p ^ 2 and therefore the anisotropy parameter 
varies according to 



0=1- 



2 + eg 
4 \ 1 + c 2 + eg 



(22) 



From its construction, the hypervirial relation is automatically 
satisfied. Because the potential is everywhere finite, the hyper- 
virial relation also implies that the velocity dispersions are ev- 
erywhere finite. In particular, assuming c + 0, the central ve- 
locity dispersions are {v 2 } = l/(p + 2) and (v£) = p/[2(p + 2)]. 
In addition, 



dp _ p 2 c(l-c 2 )r p/2 - 1 (r p - 1) 
dr ~ 8(c + r p / 2 ) 2 (l + cr p l 2 ) 2 ' 



(23) 



and therefore, provided that < c < 1, the anisotropy parame- 
ter p decreases from /3 = l-p/4atr = OtojS = 1 -p/[2(l +c)] 
at r = 1 and then increases back to P — > 1 - p/4 as r — > oo. 
In other words, in this model, the velocity dispersions are more 
radially anisotropic near the centre and the outskirts whereas 
they are relatively less radially anisotropic around the region 
of the radial scalelength. However, we also note that the model 
as a whole always possesses a more radially biased velocity 
dispersion than the isotropic model unless p > 2(1 + c). 

It has been suggested that dark matter haloes achieve an al- 
most isotropic state near the centre and become more and more 
radially anisotropic in the outer p arts, at least according t o cos- 
mological N-body simulations ( Hansen & Moore 2006). Our 
DFs are radially anisotropic and therefore better suited to mod- 
elling the outer parts and envelopes of dark haloes. They are 
unsuitable for the class of problems in which central anisotropy 
plays a critical role. 

4. An application: the distribution of y-ray sources 

iDiemand. Moore. & Stadell d2005l) have presented evidence 
from numerical simulations that dark matter may be clumped 



into mini-haloes of Earth mass and larger. They estimate that 
~50% of the total mass of the dark matter halo is bound to 
dark matter substructures. These objects may be detectable by 
virtue of the y-rays from neutra lino a n nihila tion in the very 
centres of the clumps. IDiemand et alJ J2005b also point out 
that the nearest mini-haloes will be amongst the very bright- 
est sources from neutralino annihilation and may be found ei- 
ther with the forthcoming GLAST satellite or next-generation 
atmospheric Cerenkov telescopes as high proper motion, dis- 
crete y-ray sources. 

If this idea is correct, then there are some immediate con- 
sequences. First, because of the offset of the Sun's location 
from the centre of the dark halo, the distribution of such y-ray 
sources is anisotropic and the magnitude of the anisotropy is an 
indicator of the cusp slope at the Galactic Centre. This effect is 
already well-known in studi es of halo origin of the u ltra-high 
energy cosmic rays ("e.g. jEvans. Ferrer. & SarkaJ2002l) . Let us 
use Galactic coordinates (s, I, b), where s = |s| is heliocentric 
distance, and (£, b) are Galactic longitude and latitude. Then, 
the relative number density of y-ray sources is 



F({, b)oc s z ds w(s) p[R Q + s ( ,b] 



(24) 



where Rq is the Galactocentric solar position and w(s) is the 
selection function. For bright sources, the selection function is 
proportional to the relative luminosity and so w(s) oc s~ 2 . The 
density profile is truncated at a radius at which the dark mat- 
ter annihilation rate matches the collapse timescale of the cusp. 
With this assump tion, a tiny constant density core is created 
(e.g., Tvlerl l2002l) and convergence of the integral J24I is guar- 
anteed. The overall normalization of the integral depends on 
the fraction of dark matter bound in mini-haloes, as opposed 
to smoothly distributed dark matter. The anisotropy effect is 
clearly illustrated in Fig. [5] which shows Hamer-Aitoff projec- 
tions in equatorial coordinates for two of the halo models. The 
first is a model with p oc r~' at the centre (p,c) = (2, 1/2), 
and the second is a model with p oc r^ 1 ^ 4 at the centre (p, c) = 
(1/2, 1/2). The more highly cusped the model, the greater the 
anisotropy. If there are enough detections, the magnitude of the 
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Fig. 4. The distribution of line of sight velocities and proper motions of dark matter particles towards the Galactic Centre (upper 
panels), the anti-Centre (middle panels) and the Galactic North Pole (lower panels). The left-hand panels refer to the model 
with (p,c) = (2, 1/2), while the right-hand panels refer to the model with (p,c) = (1/2, 1/2). The velocity is given in units of 
v = (GM/r ) 1/2 . 



anisotro py can be quantified by harmonic analysis (see e.g., 
lEvans et all20o3l . 

Second, if the high proper motion sources can be identi- 
fied, this will provide the first direct evidence on the velocity 
distribution of the dark matter. It is therefore useful to have 
predictions of the velocity distributions. The proper motion 
and radial velocity distribution of dark matter clumps can be 
calculated easily for our models, because the DFs are simple. 
The proper motion and radial velocity distributions in direction 
d 2 Q = cosfrdfdb is 



dv { dv h ' 



(v e , v h , I, b) = dfo, Jj dv s ds f(E, L) s 1 w(s), 



dN, 



obs , 



dv. 5 



(v s , i, b) = dfel jjj dve dv b ds f(E, L) s 1 w(s). 



(25) 



(26) 



Here, the velocity of the dark matter particles has been re- 
solved into components (y s , V(, Vj) based on Galactic coordi- 
nates while f(E, L) is the DF of the dark matter. For models 
with DFs (0, the above triple integrals are easily computed via 
Gaussian quadrature. In fact, the second moments of these dis- 
tributions are entirely analytic, as sketched out in AppendixlAl 
Fig- E] shows the proper motion and line of sight velocity dis- 
tributions for the dark matter clumps in three directions. When 



looking towards the Galactic Centre or anti-Centre in spheri- 
cal models, then the distributions of the two tangential velocity 
components vt and Vb are the same, though this is not necessar- 
ily the case in other directions. All three velocity distributions 
show very significant deviations from Gaussianity. The sources 
with highest tangential velocity occur towards the the Galactic 
Centre direction and it is therefore in this direction that GLAST 
or VERITAS may find it easiest to detect them unambiguously. 

Although these facts have been established in the context 
of our simple family of models, we argue that they are likely to 
be generic. The anisotropy effect in source positions is a conse- 
quence of the Sun's offset from the centre. The fastest moving 
sources are likely to be found in the direction in which the grav- 
itational potential well is deepest. 

5. Conclusions 

We have presented a simple family of halo models useful for 
the study of dark matter haloes. They have a simple potential 
and distribution function (DF), being composed of two terms 
of the form of £/~ 2 £ , 0+ 1 )/ 2 _ where L is the total angular 
momentum, E the binding energy and p is some constant. As 
an application, we have presented the properties of discrete y- 
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ray sources arising from neutralino annihilation in dark matter 
mini-haloes. This suggestion for the compositi on of the dark 
matter is derived the numerical simulations of Die mand et alJ 
(2005j). If true, then nearby, high proper motion discrete y-rays 
sources may be detectable by forthcoming missions such as 
GLAST. We have shown that there is an anisotropy in the po- 
sitions of such y-ray sources because of the offset of the Sun 
from the Galactic Centre. We have also provided distributions 
of proper motions and line of sight velocities of such sources. 
The best direction in which to look is towards the Galactic 
Centre {(, b) = (0, 0). There are more sources in this direction 
and they have the highest proper motion. 

This pape r also con tinues the study of hypervirial ity be- 
gun bylEdd ington and further developed bv lEvans & Anl 
(2005). In spherical symmetry, any model which has a DF that 
is a sum of terms like L P ' 2 E^ P+X ^ 2 is hypervirial. This leads to 
the derivation of a general differential equation that any spher- 
ical hypervirial system must satisfy. We have solved it to find 
the most general models, for which the potential can be written 
down in terms of elementary functions. 

We briefly note that, in the case of axisymmetry, any model 
which has a DF that is a sum of terms like L P ~ 2 E {3p+l) l 2 - where 
L z is the angular momentum component about the axis of sym- 
metry - is also hypervirial (see AppendixlBlfor more details). 
It is possible to derive the equation for the corresponding po- 
tential and demonstrate that at least two analytical solutions 
exist, corresponding to the flattened Plummer model studied 
bv lLvnden-Belll (h%2|) and a particular case among the prolate 
galaxy models o f lLakeUl98lh . In fact, only the Plummer model 
- the sole isotropic, hypervirial model - can be generalized to 
give an axisymmetric hypervirial system. It is an open question 
whether any further hypervirial generalizations of the Plummer 
model can be found. 
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Appendix A: Second moments of the velocity distributions 

It is straightforward to compute the six independent components of the velocity dispersion tensor in a Galactic coordinate system. 
These formulae do not appear to have been given before, and so we quote them here: 

{v 2 ) = (v 2 .) cos 2 b cos 2 £ + (v 2 ) cos 2 b s'm 2 £ + (v 2 ) sin 2 /? + (v x v y ) cos 2 b sin 2£ - (v y v z ) sin 2b sin £ - (v x v z ) sin 2b cos €, (A. 1) 

(v 2 ,) = (v 2 ) sin 2 /? cos 2 £ + (v 2 ) sin 2 b sin 2 £ + (v 2 z ) cos 2 b + (v x v y ) sin 2 b sin 21 + (v y v z ) sin 2b cos £ + (v x v z ) sin 2b sin £, (A.2) 

(v 2 ) = <v 2 )sin 2 ^ + <v 2 )cos 2 ^-<v x v v )sin2^ (A.3) 

9 i 9 9 

(v s Vb) = —(v x ) cos b sin b cos £ - (vo cos b sin b sin £ + (vp cos b sin b - (v x v y ) sin 2b cos £ sin £ 

— (v x v z ) cos 2b cos £ - (v y v z ) cos 2b sin I , (A.4) 

{v[V s ) = (v 2 - v 2 ) cos sin £ cos £ + (v x v y ) cos b cos 2€ + {v x v z ) sin Z? sin ^ - (v y v z ) sin b cos (A.5) 

(vtVb) - {v 2 - v 2 ) sin b sin /" cos £ - (v x v y ) cos cos 2£ + {v x v z ) cos b sin £ - (v y v z ) cos b cos £ (A.6) 

This gives the observables in terms of the second moments referred to Cartesian coordinates. The latter can be related to the 
second moments in spherical polars (in which the tensor diagonalizes) via: 

/ 2\ X 2 i 2\ , y 2 + ^ / 2\ / 2\ y 2 , 2\ , X 2 + Z 2 , 2\ i 2\ ^ / 2\ , x2 + 3^ , 2\ / a t\ 

<vj> = -2<v r ) + r2 <V g ), <v y ) = ~^{v r ) + r2 {vp, {Vp = -2<v r ) + — {vp, (A.7) 



{v x v y ) = ^(<v 2 )-<v^>), <v,v z > = ^(<v 2 >-<v 2 >), < VjVz > = Z£(< v 2 >-<v 2 >). (A.8) 



Here, {vp = (v?) = (Vj)/2. Hence, by substituting in the known velocity dispersions in spherical polars (Eas.l20land l2U . we 
obtain the second moments of the observable distributions (the line of sight velocities and the proper motions). 

Appendix B: Axisymmetric hypervirial potentials 

Upon the discovery of the new hypervirial potential, one may ask the question: does there exist a simple extension of the spherical 
hypervirial potentials into axisymmetry? A possible starting point of the investigation is the examination of axisymmetric DFs 
that are obtained by replacing the L-dependence of a spherically symmetric DF with an independence. Let us suppose that the 
DF of an axisymmetric system is given by 

f(E, L z ) = C\L z \ 2 "E' n - 3/2 , (B.l) 

where L z = Rv^ is the angular momentum along the symm etry axis, and n > —1/2 and m > 1/2. Then, a simple integration leads 
to the expressions for the corresponding density (c.f., Fricke 1952; Deionahe 1986; Evans 1994) 

= ^3/2^^ + 1/2)^-1/2)^ 

F(n + m + 1) 
and the second-order velocity moments 

p{v 2 R ) = r^c r(n+1/2Mm - 1/2) R 2 = ^ ; 

Y(n + m + 2) n + m+ 1 

p<v 2 > = 2" +5 / 2 J rC r(W+ r f /2)r(m " 1/2) ^y +m+1 = (2n + l)p<v 2 >. (B.3) 
v Y(n + m + 2) 

Note that p{v 2 R ) = p(v 2 ) = p(v 2 ) = p(Vg) since the DF is symmetric with respect to vr <-> v, exchange. Finally, we find that, for 
the system described by the DF of Eq. (IB. li . 

P<v 2 ) = P (<v|> + <v 2 > + <v 2 >) = (2n + 3)p<v 2 ) = p^, (B.4) 

and thus that the steady-state virial theorem holds without the boundary terms only if m — 3n + 5, for which the system becomes 
hypervirial. Here, the density of the system (Eq. IB.2l is axisymmetric, but there is no net rotation ({vp = 0), as the DF ( IB. 1 1 is 
symmetric with respect to L, «-> —L z exchange. To rectify this, we can amend the DF JB. 1 i by 

f(E, L z ) = C [1 + f(Lj)] |L : | 2 "£'"- 3/2 (B.5) 

where £(L C ) is an arbitrar y odd function of L 7 fi .e, £(-L z ) = -£(L,)] that is bounded by - 1 < £(L Z ) < 1 . Note that, on the ground 
of statistical mechanics, Deionghe ( 1986l ll987l) advocated the form of g(L z ) being tanh(^L-) where k is a parameter that depends 
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on the total angular momentum along the symmetry axis. Physically, this action of adding a certain odd function to the even 
distribution function corresponds to switching the direction of rotations for a certain fraction of stars of given \L Z \ and E. It is 
straightforward to show that, for the DF of Eq. ( IB.5I . Eqs. \B.2\ and iB.'H are still true while 

, I^-m . .../ r 2 \'"- 1/z 



^ 2 >" + '" +1/2 f^fylftj f (1 - t) m - ljl df 



2 n+i xC 
2m — 1 



< 2 n + 2 c r (» + l)r(w - 1/2) 2 +m+1/2 _ 1/2 r(n + l)r(n + m + 1) 1/2 
r(n + m + 3/2) V Y(n + l/2)T(n + m + 3/2) P ^ ' 

In other words, one can choose %{L Z ) to meet the constraint on the net rotation (v^). Here, the inequality becomes the equality 
if %(L Z ) = 1 for L z > 0, which corresponds to the situation when every star rotates in the same direction. For n > -1/2 and 
m — 3n + 5, it is straightforward to established that (v^) 2 < iff/2. 

Systems given by the DF JB.lt are somewhat unrealistic since the density along the symmetry axis is either zero (« > 0) or 
infinite (n < 0) except the n — case that is just the isotropic (and spherically symmetric) Plummer model. However, analogous 
to the procedure in Sect. |2 we may obtain more realistic systems by summing DFs of the form of Eq. JB. It , one of which 
corresponds to that of the isotropic Plummer model, that is, 

f(E, L z ) = C E 1/2 + ^ d\L z ^E^ in + f odd (E, L z ) (B.6) 

i 

where each «, should be positive in order to avoid the divergent behaviour of the density along the symmetry axis. Here, f dd(E, L z ) 
is an arbitrary function that satisfies the condition that f dd(E, -L z ) = —f d&(E, L z ) and chosen so that the f(E, L z ) is non-negative 
for all accessible E and L z . It is a simple exercise to find expressions for the density 

P = + V A* 2 V- 5 ; D = ^C ; D, = T'^C ^ + ^f^ + ^ . (B.7) 
t-r* 2 n/z r(4«/ + 6) 

and the second-order velocity moments 

p(v 2 } = £y + y — ; P (vi) = £y + y (2 "' + l) °' i^y 4 ** (B .8) 

^ s 6 4-* 2(2n,- + 3) ^ 6 2(2n,- + 3) 

It is also easy to establish that the system is in fact hypervirial, 

p<v 2 ) = p (2<v 2 ) + <v 2 )) = ^ + Z| ^V*" = Pf • ( B - 9 ) 

i 

However, finding the corresponding potential is actually a non-trivial problem as it is tantamount to solving Poisson's equa- 
tion, which now becomes a second-order partial differential equation: 

LI + UnO^] = -4jzD ^ 5 -4k y D ! r 2 " i wfrotf** 5 (B.10) 



r 2 dr I <9r / r 2 sin # dr \ 36* , 
in the spherical polar coordinate system or 

ia(*ffi) + £-^-<«?***~ 

in the cylindrical polar coordinate system. Even after factoring out the scalefree contribution (c.f. lToomrell982l) . 



a? + cosh ^ " 4 = - A * D °* - 4 "2 fl '^ < B - 12) 

where <^> = r^ 2 i/f, f = In r, and cos 6* = tanh f or 

where </> = R l ^ 2 i//, these are not obviously separable in either coordinate system. Instead, here, we take a rather different approach. 
Inspired by Eq. H51 . we first consider the axisymmetric potential of the form of 

>ff = l -r, rr (B.14) 

r [1 +c(ff)rPl 2 + rP] 1 'P 
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where c(9) is some function of the latitude 9. Then, the corresponding density can be found to be (G = 1) 



p = __Lvfy = (p + i) 

4jt 



1 - 



1 dc 
pd6» 



r P-2^2 P+ i + I 
P 



Id/ dc\ p(o + 2) 
sin0— + — -c 



sin6»d6>\ 60 



r p/2 - 2 i/r p+1 . (B.15) 



Now, we find that one of the necessary condition for the system to be consistent with the hypervirial DF is that c(9) is the solution 
of the inhomogeneous Legendre equation 



Id/ dc 
sin 



sin 9 6.6 



where C is a non-negative constant. A simple trial reveals that en = 4Cp _I sin'^ 2 6> is its particular solution, and so the general 
solution is 

AC 

c = AP p/2 (cos 9) + BQ p/2 (cos 9) + — sin p/2 (B. 17) 

where Pi(x) is the Legendre-P function (polynomial if I is an integer), Qi{x) is the Legendre-Q function, and A and B are the 
integration constants. Finally, the system can be hypervirial only if 

1 - — - f-— V ccshrt (B.18) 

where c is given by Eq. flJ3. 17b . It is rather obvious that we only need to examine the cases when p/2 is an integer and B = 0. 
After a further examination, we find that, only possible cases are p — 4 

c(9) = 3 cos 2 9 - 1 + C sin 2 9 = 2 + (C - 3) sin 2 0; ^ = - - — - -g i - _ (B.19) 
and p = 2 

c(9)= A cos 9 + 2Csin9; # = 1 (B.20) 

[1 + r 2 + Az + 2CR] 1 ' 2 

where A 2 < 4(1 -C 2 ). The corresponding densities for both cases reduce the form of Eq. iB.H . In fact, the potential corresponding 
to both solutions are already known! 

The first (p = 4) is the flattened Plummer model of lLvnden-Belil ill 9621) . 

<foa _ <foO - 2e) / <A \ 5 5^)6(2 - e) [R\ 2 I # \ 9 

V [(r 2 + a 2f_ 2b 2 R 2 }^ ' 9 xGa 2 Uo/ xGa 2 \a) \iff ) 

where < e = b 2 /a 2 < 3/2. Comparing this to Eq. ( IB.7L we find that the DF 

A-flW-W^ + G^ «-%'<£,- c ' = ^fvT (B - 22) 

can build the potential-density pair of the flattened Plummer model. The second-order velocity moments are 

= <A ( 2 (3-2 C ) / 1\ 6 + <A 2 e(2-e) /^ 2 /^\ 10 _ = ^(3-2.) / ^ + 3^(2 - f ) - vio ^ ^ 

A 6jiGa 2 \iAo/ 2nGa 2 \a) \tl/o) ' * 6jrGa 2 V^o/ 2jiGa 2 \a/ \tAo/ 

and the system is hypervirial, as already noted bv lLvnden-Be1illll962l) . 

The second (p — 2) is a particular case of the generalized Plummer model devised bv lEakel lll98ll). for which 

■ *!>oa = 3iAo(l - c 2 ) / iff \ 5 /«\/'A\ 3 rR?4 ^ 

^ [r 2 +a 2 +2M] 1 / 2 ' P 4jxGa 2 Uo/ 4jrGa 2 U/U<J ' 

where < e = Z?/a < 1. This reduces to Eq. ( IB. 191 after an arbitrary translation along the symmetry axis {R = 0). In fact, this is 
singular on the entire R = axis and so best represents prolate galaxies. The DF is 

f evw (E,L z ) = C E 1 ' 2 + C 1 S(L z )E 2 ; C = ; d = ^-^-j. (B.25) 

7jt j Ga 2 ^ 8jt 2 Gai^o 

Note that the Dirac-5 distribution is the limiting case of Eq. J B . 1 i as in lim„^(_i/2)[|L,| 2 "/F(« + 1 /2)] = 6{L Z ). The second-order 
velocity moments are 

„ 2n _ ^(1-e 2 ) /<aV , Ke ^j/£j 4 . „ / „ 2x _ ^d- e 2 ) / l A 



^> = -W- + i6^UJfe) ; ^ ) = i^fei ' (B - 26) 
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and so this model is hypervirial. 

The existence of the flattened Plummer models poses the question as to whether there are axisymmetric generalizations with 
simple DFs for all the hypervirial models. We argue that this is not the case because the Plummer model is the only isotropic 
hypervirial model. Any DF of f(E, L) of a spherically symmetric system has degenerate velocity dispersions in the tangential 
plane, while any DF f(E, L z ) of a axisymmetric system has the velocity dispersions within the meridional plane degenerate. 
Clearly, only the isotropic case can have both the tangential plane and meridional plane as planes of symmetry of the local 
velocity ellipsoid. This, therefore, implies that, if the flattening method, however it is chosen, is gradual in the sense that it 
includes the spherically symmetric case as a particular case, the corresponding DF should necessarily reduce to not only a 
spherically symmetric but also an isotropic DF of f(E). In other words, only those spherically symmetric potentials with a 
simple isotropic DF can have an axisymmetric family of its generalization with continuous values of flattening parameters. This, 
of course, does not rule out the possibility that with a very special values of flattening parameters, the corresponding DF may be 
simple but there would not be any continuous transitions to the spherically symmetric case while maintaining the simplicity of 
the DF. 



Appendix C: General solutions of generalized Lane-Emden equation 

The differential equation in Eq. w hen the sum contains only a single term is a particular case of a family of differential 
equations investigated by Goenne r & Havasl J2QO0l) . which they called the generalized Lane-Emden equation. I n fact, Eq. (|6j 
corresponds to the special case of parameter combinations noted by them (see eq. 1 1 of lGoenner & Havaslfc oOO) that permits a 
ra tional transformation of v ariables that leaves the differential equation form-invariant. The one -parameter so lution family (eq. 18 
of lGoenner & Havasl2000l) is the direct generalization of the Schuster-Em den integral (see e.g . lHoredil986l) and corresponds to 
the solution for the case of A — in Eq. (1101 . which leads to the models of E vans & Anl d2005l) . 
If we consider a transformation, 77 = ip s , Eq. fllOi with a single term in the sum further reduces to 

= AsW~ 2/s + jl 2 ~ BsW +2p/s - (CD 

We note that the differential equation of the form of (y') 2 = f(y) where f(y) is a cubic or quartic polynomial of y can be solved 
with standard elliptic functions. Here, for AB + 0, the right-hand side of Eq. dC.U can become a cubic of quartic polynomial of 
r\ if s — ±1 (with p — 1 or p — 2) or s — +2 (with p = 1/2 or p = 1). In other words, the two-parameter general solutions of 
Eq. or equivalently the solutions of Eq. dlOt with A + can be wri tten down in closed form using elliptic functions if the 
sum contains a single term with p =1/2, 1, or 2 (G oenner & Havasl2000l) . 



C.1. Solutions in terms of Weierstrass-P functions 

Let us be reminded that the Weierstrass-P function is the canonical elliptic function (i.e., complex bi-periodic meromorphic 
function) with second-order poles. It is usually defined in terms of the sum over the lattice points in complex plane with two 
half-periods. However, for our purpose, it is useful to consider the differential equation 



dz 



P(z;g2,g3) 



4[p(z; gz, gi)] - g2p(z; g2, gl) -g3 



(C.2) 



that is satisfied by p(z',g2,g3)- Here, g2 and #3 are usually referred to as elliptic invariants. We further note that the differential 
equation dC.2t does not explicitly involve the independent variable so that one integration constant is given by an arbitrary 
translation, that is, if ip(t) is its solution, then <p(t + c) where c is an arbitrary (complex) constant is also the solution. However, 
Eq. dC.2t is a first order differential equation, and thus, p(z + c;g2, gi) is indeed the complete specification of its general solution. 
In addition, the above differential equation also indicates that p-function is homogeneous, i.e., 



p{Az\gi,g3) = A 



l(z;A 4 g 2 ,A 6 g 3 ). 



(C3) 



The right-hand side of Eq. dC.H becomes a cubic polynomial if i) B + 0, s — 2p (p — 1/2 or 1) or ii) A + 0, s — -2 (p — 1 
or 2). For these cases, linear transformations 77 = A(£ + r) with prop erly chosen A and j can reduce Eq. dC.U to the form of 
Eq. dC.2l >. After a bit of algebra, we find solutions of Eq. JlOi (see also lGoenner & Havasl200(il) : 



1 



2 2 -3 



-2 2 p f + c; 



1 



2 6 -3' 



AB 2 



1 



2 9 .3: 



for B + 0, p = 1/2, 



(C4) 



* =B 



2 2 -3 



, 1 AB 

t + c; 2 2 AB + -= — , 

2 2 -3 3 



2 3 -3 3 



for B + 0, p = 1, 



(C5) 



-2 1 
-2 1 



f + c; 2 2 AB + 



p f + c; 



2 2 -3 
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1 AB 1 \ 1 

for A #0, p = l, 
for A + 0, p = 2. 



2 2 -3' 3 2 3 -3 3 / 2 2 -3 
1 \ 1 



2 2 A 2 fi 



2 3 • 3 3 / 2 2 ■ 3 
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(C6) 
(C7) 



C.2. Solutions for p = 1 in terms of Jacobi elliptic functions 

The Jacobi elliptic functions are defined through the inverse function of the elliptic integral of the first kind, i.e., 

sn(u, k) — sin $ ; cn(u,k) — cos</>; dn(u,k)- -yl - k 2 [sn(u, k)] 2 

where 

dt f sn(ia) dv 



m = F(0 



r dr f 
J o VI - £ 2 sin 2 f J » 



(1 _ v 2)l/2(! _ £2 v 2)l/2 



•/t 2 sin 2 f j!i 

and is the elliptic modulus. Remaining elliptic functions are defined as 



ns u = sn u; 
sc u = sn m nc m; 
cs m = ns u cnu: 



nc m = cnu; 
cdu = cnu ndu; 
dc u = ncu dnu: 



ndu = dnu; 
dsMEdnw nst<; 
sdM = ndM snM. 



(C8) 



(C9) 



(CIO) 



(Cll) 



Here, the common elliptic modulus k is omitted following usual practices. The derivative of sn u is given by 
d 

— sn u — cn u dn u 
du 

and the derivatives of remaining elliptic functions can be found using their respective definitions and Eq. jC.ll> . Then, it is 
straightforward to show that each Jacobi elliptic function is a particular solution of second-order differential equations such that 

(sn u)" - -(1 + k 2 ) sn w + 2A: 2 sn 3 «; (cn u)" = (2k 2 - 1) cn u - 2k 2 cn 3 u; (dni<)" = (2 - k 2 ) dnu — 2dn 3 w; 

(ns u)" = -( 1 + k 2 ) ns u + 2 ns 3 u; (ncu)" = (2k 2 - 1 ) nc u + 2( 1 - k 2 ) nc 3 u; (nd u)" = (2 - k 2 ) nd u - 2( 1 - k 2 ) nd 3 u; 

(2k 2 - l)sdu-2k 2 (l -k 2 )sd 3 u; (scu)" = (2 - k 2 ) sc u + 2(1 +k 2 )sc 3 u; ' 



(cd u)" = -(1 + k 2 ) cd u + 2k 2 cd 3 M ; (sd u)' 
(dcu)" = -(1 +/fc 2 )dc M + 2dc 3 M ; (dsu)" = (2k 2 - l)dsn + 2ds 3 w; 



(cs u)" — (2 - k ) cs u + 2 cs u 



where the primed symbols indicate the differentiation with respect to the argument u and the simplified notations for power (e.g.. 
sn 3 w = [sn(u,k)] 3 etc.) are used. 

For Eq. (|9j when the sum contains a single term with p = 1; 



dV 
dt 2 



t 
4 



2B<p 3 



(C.13) 



where B = 2jtD, we first note that it does not involves the independent variable explicitly, and therefore that, if <p(f) is its solution, 
then (f(t + c) where c is an arbitrary constant is also the solution. That is, one of the integration constants of the general solution 
of Eq. ( IC. 1 31 is related to arbitrary translation of a particular solution. Next, let us think of linear transformations of variables, 
(f — Aip and t - a(t — to) where to is an arbitrary integration constant while A and a are constants to be specified. With these new 
variables, Eq. (IC. 13i is transformed to 

d 2 (p 



_L~_ 2 ^~ 3 

dt 2 Aa 2>P a 2 V ' 



(C14) 



which can reduce to any of the differential equations in Eq. dC.12t by proper choices of the constants. Note that there are two 
constants, A and a, to be specified but Eq. AC 12b contains only one parameter k 2 . This leaves one arbitrary degree of freedom, 



which essentially provides with the remaining constant of integration. Then, we find the general solutions of Eq. (IC . 1 3 i 



1 +c 
~~8B~ 

8A 



1/2 



cn[a + (t-t ),k+] ; 



4A\ 1/2 

— I sd[a + (t-t ),k + ] ; 



1+c 



1/2 



nc[a+(t- t ),k+] ; 



=4~- 

\ AB 



M 2 



ds [a+(t- t ),k+] ; 



1 +c 



1/2 



dn [a-(t - t ), k-] ; 



8A 
1+c 



1/2 



nd[a-(t-t ),k-] ; 



(C.15) 



(C.16) 



(C.17) 
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(8 A \ il~t~c\ 

— I sc[a-(t-t ),k-] ; <p = ±[— g^l cs [»_('- to),*-] (C.18) 

where 

c = {l+64AB)^>0; al = °-, al = I±£ ; ^ = 1 |j + Ij = ^ (C19) 

and A and fo are two integration constants. In particular, A is chosen to correspond to the constant in Eq. (11 Ot . Here, all solutions 
are equivalent to one another in a sense that one can be transformed to another with a proper choice of the integration constant 
(which is, in general, allowed to be complex). 

We note that, in principle, the solutions for p = 1/2 and p — 2 can also be written in terms of Jacobi elliptic functions. 
However, typical reduction procedures involve more complicated transformations as well as the determination of constants though 
the solutions of cubic or quartic polynomials. 



C.3. Solutions expressible using elementary functions 



The p-function reduces to elementary functions if the cubic equation f(z) = 4z - giz - gi — has a degenerate solution or 
equivalently g\ = 21 g 2 . If the degenerate solution is zero, then g2 = gj = and p(z; 0, 0) = z~ 2 . If e + is the nonzero degenerate 
solution of the cubic equation, since the quadratic coefficient of the cubic polynomial f(z) is null, the other solution of the cubic 
equation is -2e while g2 = I2e 2 and gi 

p(z;12e 2 ,-8e 3 ) = 



-8<? . Then, we find that 



3ecoth 2 [(3e) 1/2 z] -2e = 3ecsch 2 [(3e) 1/2 z] + e 
(-3e)cot 2 [(-3e) 1/2 z] + (-2e) = (-3e) csc 2 [(-3e) 1/2 z] - (-<?) 



if e > 
if e < 



Similarly, Jacobi elliptic functions reduce to elementary functions if k — or k — 1 . For example, 

sn(M, 0) = sin u; cn(u, 0) = cos u; dn(M,0)=l; 
m(u, 1) = tanhM; cn(M, 1) = secht<; dn(M, 1) = sechw. 



(C.20) 



(C21) 



Using this, it is possible to find some nontrivial solutions of Eq. (|6j which can be written using elementary functions for the 
case when the sum contains a single term with p = 1/2, 1, or 2. 1) For p = 1/2, they are 



1 



6Br 1 / 2 ' 



1 



< 2 ( r l/2 + r l/2)2 



ro 



Br, 



2 ' 



<A = - 



1 



UBr 1 ' 2 



1 + 3 tan 2 



-lnl- 

ro 



where 3B = 8jiZ) + 0. II) For p = 1, 

1 1 r„ 



0r = ± 
lff = ± 



(8Br) 1 / 2 
1 



-A 



-(Br y/ 2 r + r 



■ tan 



1 



V8 



;ln 



(-%B) l l 2 r l l 2 
where B = 2nD. Ill) For p — 2, 
1 



4r 



for B > 0, 
1 



'0 



" \-Byi 2 r]j 2 r - r 



<fr = 



"(12B)'/V/ 2 ' 



1 



\-ByiV' 2 



(r 2 -rD 



1/2 ' 



(3fi)l/4 r l/2 



1 + 3 tan 2 



-1/2 



for B < 0, 



^ = ± 



'o 



+ r 2 ) 



1/2 



(-B)l/^ 2 (_ r 2 + r 2) 



1/2 



for B < 0, 



(C.22) 

(C.23) 
(C.24) 

for B > 0, (C.25) 
(C.26) 



where 3B = AnD. For all, ro > is an arbitrary constant of integration. In addition, if B — 0, there are common solutions for all 
b 



i[s = a + 



for B = 



(C.27) 



where a and are arbitrary constants. These solutions include the scale-free so lution (i/r <x r ~ 1 ^ 2 w ith the proportionality co- 
efficient uniquely determined by p) and the generalize d Plummer potentials of IE vans & Anl J2005). With proper scaling, the 
oscillating solution of p — 2 case reduces the ISrivastaval \ 1 9621) solution of Lane-Emden equation of index 5 in three dimension. 
However, except for the generalized Plummer potential, none of these is physical. 



